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Abstract. In this paper we investigate the critical collapse of an ultrarelativistic 
perfect fluid with the equation of state P = (T — l)p in the limit of T — > 1. We 
calculate the limiting continuously self similar (CSS) solution and the limiting 
scaling exponent by exploiting self-similarity of the solution. We also solve the 
complete set of equations governing the gravitational collapse numerically for 
(r-1) = 10- 2 ,...,10" 6 and compare them with the CSS solutions. We also 
investigate the supercritical regime and discuss the hypothesis of naked singularity 
formation in a generic gravitational collapse. The numerical calculations make use 
of advanced methods such as high resolution shock capturing evolution scheme 
for the matter evolution, adaptive mesh refinement, and quadruple precision 
arithmetic. The treatment of vacuum is also non standard. We were able to 
tune the critical parameter up to 30 significant digits and to calculate the scaling 
exponents accurately. The numerical results agree very well with those calculated 
using the CSS ansatz. The analysis of the collapse in the supercritical regime 
supports the hypothesis of the existence of naked singularities formed during a 
generic gravitational collapse. 



PACS numbers: 04.20.Dw,04.25.Dm,04.40.Nr,04.70.Bw,02.60.-x,02.60.Cb 

1. Introduction 

Since the discovery of critical phenomena in gravitational collapse of a massless scalar 
field by Choptuik in 1993 [I] researchers have observed critical phenomena in many 
different matter models. In this paper we focus on perfect fluid, in particular the 
ultrarelativistic perfect fluid with the equation of state 

P = (T-l)p. (1.1) 

The study of critical phenomena in ultrarelativistic perfect fluids has quite interesting 
history. It started with Evans and Coleman's calculations 2 who were able to obtain 
critical solution for the radiation fluid (T = 4/3) by employing the CSS ansatz and 
also by directly solving the PDE's. Later, Koike et al [Hj used perturbation analysis to 
calculate the scaling exponent of the radiation fluid collapse. Maison ^ generalized 
the calculations for 1.01 < T < 1.888 and similar analysis was carried out by Hara, 
Koike and Adachi |5J. Both speculated that the solutions might not exist for T > 1.89. 
Fully numerical treatment of the collapse for 1.05 < T < 1.5 was performed by Evans 
and Perkins UJ. The inability to obtain the solutions for T > 1.89 lead to various 
ideas about what happens in that regime. Some proposed that the solution changes 
from CSS into DSS (discretely self similar) or a mixture of both. It was also proposed 
that the critical solution might be of type I and not type II. 
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The confusion was definitely resolved by Neilsen and Choptuik 7 who 
demonstrated by directly solving the PDE's that nothing special happens for T > 1.89. 
They carried out the calculations all the way up to T = 2. Moreover, they were able 
to obtain the solutions for T > 1.89 with the CSS ansatz. The key to success was to 
use arbitrary precision arithmetic. 

The purpose of this paper is several. We confirm the results of the perturbation 
calculations with the CSS ansatz from previous works for values of T close to 1 
and, more importantly, we confirm by direct numerical collapse calculations that the 
perturbative calculations give the correct solutions for that regime. We obtain the 
exact numerical value of the limiting scaling exponent that has never been calculated 
before. 

One of the reasons for studying the T — ► 1 limit is to verify the hypothesis made 
by Harada and Maeda that for T < 1.0105 the final state of a supercritical collapse is 
not a black hole but a naked singularity [5] . It is known that the critical solution 
contains a naked singularity but this is only achieved by fine tuning of the critical 
parameter. The naked singularity formed during the supercritical collapse would not 
require any fine tuning and therefore is much more interesting. 

Critical collapse possesses a unique set of challenges for numerical treatment. 
During the evolution the dynamics occurs on ever decreasing length scales and the 
magnitudes of the fluid state variables increase by many orders. Therefore, advanced 
numerical techniques have to be used to perform these calculations. 

For this purpose we develop a general adaptive mesh refinement (AMR) algorithm 
which is capable of performing the critical collapse calculations. This algorithm is 
general and does not rely on prior knowledge of the dynamics. We also use innovative 
treatment of vacuum that does not use a low density artificial atmosphere — a method 
widely used to ameliorate numerical instabilities at the matter-vacuum boundary. 

The use of quadruple precision arithmetic and the above-described advanced 
numerical techniques allow us to obtain the solutions and scaling exponents with 
higher accuracy than in previous works. 

The organization of the paper is as follows. Section [21 reviews the basic equations 
of general relativistic hydrodynamics in spherical symmetry. In section [3] we sketch 
the ideas behind the CSS ansatz and present the equations and their solutions for 
the limiting case r — > 1. Similarly, in section^ we briefly review the perturbation 
theory and present the equations and their solutions for the relevant eigenmodes in 
the limiting case. In the later two sections we closely follow the approach described 
in ^U] where all of the missing details can be found. Section describes the 
method of obtaining the critical scaling exponent from numerical calculations and a 
way of estimating the errors. It summarizes and compares results obtained from both 
approaches. We also compare our numerical results from supercritical calculations 
with the expected CSS solutions. In section|H|we describe in more details some of the 
features of the numerical algorithms such as AMR and vacuum treatment. Sectional 
summarizes the results and concludes. 

Throughout this paper we use geometrized units with G = c = 1. 

2. Hydrodynamics in spherical symmetry 



In our calculations we use spherical polar coordinates in which the metric has the form 
ds 2 = -a(t, rfdt 2 + a(t, r) 2 dr 2 + r 2 dfl 2 . (2.1) 
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It is convenient to define an auxiliary variable m(t,r) 

l(i L 

2V o(i,r) 

which could be interpreted as the total mass enclosed within the radius r. In the 
vacuum regions we demand that the metric reduces to the Schwarzschild one, i.e., 

a(t,r) = -^. (2.3) 

The stress-energy tensor for perfect fluid has the form 

= {p + P)u»u u + Ptf» , (2.4) 

where P is the isotropic pressure, u M is the 4- velocity of the fluid element and p is the 
total energy density as seen by Eulerian observers (see later). The energy density p 
can be subdivided into two contributions 

p = Po(l + e), (2.5) 

with po being the rest mass density and e being the specific internal energy. In 
the ultrarelativistic limit the internal energy contribution dominates the total energy 
density, i.e., 

e>l. (2.6) 

The Eulerian observers' 4- velocity u M = n^, where n M is the unit normal vector 

to the spatial hypersurfaces defined by the constant time t slices. These observers see 
the fluid moving with the 3-velocity 

v r 

v r = — t ■ (2.7) 
It is useful to define a velocity v as 

v = av r (2.8) 
since the Lorentz factor is then a simple function of v 

W = au t = J_ . (2.9) 
V 1 — v z 

The equations of motion for the fluid can be derived from conservation laws 

T^ v . v = , (2.10) 

J% = , (2.11) 

where 

= pou^ . (2.12) 

These equations must be supplemented with an equation of state which for perfect 
fluid has the form 

P=(T-l)p e. (2.13) 

In the ultrarelativistic limit described by l|2.6|) the equation of state has the form 

P=(T-l)p. (2.14) 

The equation of state l|2.14|l is the one we use. The set of primitive variables (p, v) 
characterize the state of the fluid completely. Note that since po does not enter the 
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equation of state l|2.14[l the equation (|2.11|) (baryon number conservation) is no longer 
needed for the dynamical evolution of ultrarelativistic fluid. 

We define a complementary set of conservative fluid variables (S, E) 

S=(p + P)W 2 v, (2.15) 

E= (p + P)W 2 -P . (2.16) 
The equations of motion for the fluid have the form 

S+-^ [r 2 X(Sv + P)]' = 9, (2.17) 

E+\[r 2 XS}' = , (2.18) 



where 



X = - (2.19) 

a 



and the source term '5 has the form 

/ ttz \ Tn 2q?_P 

* = (Sv — E) ( SnaarP + aa— + aaP— H . (2.20) 

V t J r ar 

The dot denotes differentiation with respect to the coordinate time whereas prime 

denotes differentiation with respect to the radial coordinate r. 

As has been argued in |7| in numerical applications it is convenient to work with 

a derived conservative quantities 

<f> = E-S, (2.21) 

U = E + S. (2.22) 
The equations of motion for these new variables have the form 

q+l[r 2 Xf]' = S, (2.23) 

with 

* - \ |(n-$)(i-«)-p ) ' ^ -{-*)■ {2 - 2A) 

The equations for the geometry can be obtained from the non-trivial component 
of the momentum constraint 

a = -4Trraa 2 S , (2.25) 

the polar slicing condition 

^ = a 2 (4nr(Sv + P) + ™) , (2.26) 

and the Hamiltonian constraint 

a' o I rn\ , 

— = a 2 (<inrE~^) . (2.27) 

a V r A J 

In the numerical calculations we use only equations (|2.26|l and l|2.27|l . i.e., we use 

a fully constrained evolution for the geometry. 
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3. CSS solutions 

Under the assumption that the critical solutions are continuously self similar we can 
transform the system of PDE's (j2.23|) and l|2.25|) - [|2.27|) into system of ODE's that can 
then be solved numerically. Once the critical solution is found it is possible to solve 
for linear perturbation modes. The relevant mode is the mode that has the eigenvalue 
with the largest real part. There exists a simple relation between the eigenvalue of 
the relevant mode and the scaling exponent 

7=i, (3.1) 

where n is the eigenvalue of the relevant mode (in our case it is a real number). 

Let us now briefly summarize the procedure of finding the CSS solutions. The 
details can be found in [S] . We can rewrite the equations i|2.23ll , (|2.25l) (|2.27|l in terms 
of new variables adapted to the CSS symmetry 

s = -ln(-t«) , (3.2) 
x = ln(~). (3.3) 

The time coordinate t* is fixed by the requirement that the collapsing CSS solution 
reaches the origin at time i* = and that the sonic point is always located at x = 0. 
Using a new set of variables 

(3.4) 

(3.5) 
(3.6) 





r, 


N = 






ae x ' 


A = 


a 2 , 


UJ 


4-7rr 2 a 2 p 



we can write the system 



^ = i _ A + Ml±£r x > 

A 1-v 



1-A+ ^ ± ^'\~ 1 ' U , (3.7) 



AT* = _ 
N 

A, A T 2TNvu 



A-(2-r>, (3.8) 

(3.9) 



A A 1-v 2 ' 

(T-l)v (1 + Nv) T(N + v) 

* -U.s + -W.x + 5-^, x 

UJ UJ 1 — V 

2>(2 — T) 2 + r 
= -— - — -Nv —ANv + (2 - T)Nvuj , (3.10) 

(r - i)v r n + v m + Nv) 

uj s + xV . + (1 - 1) uj x H — v x 

UJ 1 — v z UJ 1 — v z 

7T — fi 9 — 3F 
= (2-T)(r- 1)Nuj-^— N+ AN . (3.11) 

These equations are not independent — the equation 1)3.9(1 is automatically satisfied 
if the rest of them are providing the solutions satisfy the boundary conditions 

A(s, -co) = 1 , 
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v(s, -co) = u/(s, -co) = . (3-12) 

For a CSS solution all of the variables (N, A, v, uj) are functions of x only therefore 
the equations (|3.7|) - (|3.11|) reduce to a set of ODE's 

A x 2uj(1 + (r - l)v 2 

~f = l - A+ i\ 2 ' . (3-13) 

TV 

_^ = -2 + A-(2-I> , (3.14) 

A iX _ 2TNvuj 
~A ~ ~ 1 - v 2 ' 



(3.15) 



+ = - + (2 - P)^ , (3.16) 

a; 1 — v 2 2 

,N + v T(l + Nv) , w , 7r-6 2-3r , 

r - 1 uj, c + — — ^-v >x = 2 - r r - 1 iv w — — n + -^—an . 3.17) 

From l|3.13|l and l|3.15|l an algebraic equation can be formed 

(1 - A)(l ~ v 2 ) + 2uj(1 + (r - l)v 2 = -2TNvuj , (3.18) 

which can be, in principle, used to eliminate one of the variables. We will label the 
CSS solutions as (N ss , A ss , v ss , w ss ). 

We are interested in solutions with T close to 1. In what follows it is convenient 
to define a new parameter 



k = Vr - 1 . (3.19) 

Our strategy is to expand the CSS solutions in powers of k. We can deduce the 
leading terms by looking at solutions for k close to zero. Our ansatz is 

N ss - , (3.20) 
k 

A ss (x) = 1 + A{x)k 2 , (3.21) 

w ss {x) = Lb{x)k 2 , (3.22) 

v ss {x) = v(x) k . (3.23) 

Substituting the expressions (|3.19() - (|3.23() into the equations (|3 - 1 3|l — (|3 . 1 8|) and keeping 
only the leading terms we obtain the equations for lu and v 

(1 + N Q e- x v)Q >x + N e- X u)v, x = , (3.24) 
Nn - - - 

-^■u x + (e x + N v)v x = -N (Cj + N Cuve- x - 2) . (3.25) 
w 

The algebraic equation (|3.18() yields an expression for A 

A = 2w(l + N e~ x v) , (3.26) 

which we used to eliminate A from the equations. 

The scaling relations (|3.20() - (|3.23|l are those of Newtonian CSS solutions. Indeed, 
the equations (|3.24|) - (|3.25|) are exactly the CSS Newtonian equations as can be seen 
by comparison with equations (3.5)-(3.7) of |llj after the transformation x = —e x /N . 

The numerical solutions of the equations H3.24fl - 13.25fl is shown in figure 2] In 
this solution the velocity crosses the zero line once and it is the critical solution. It is 
the so-called Hunter's solution of type (a) [T^] . 
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Figure 1. 


Plot of the limiting solutions of N = Noe~ x , A(x), ui{x), v(x). 



4. Perturbations of the CSS solutions 

Again, following jS] we calculate the relevant perturbations using a linear expansion 
H(x, s) = H ss (x) + e 

) , (4.1) 

where H(x,s) is one of {log(TV), \og(A), log(iy), v} and h val (x,s) has the form 

h var (x, s) = h p (x)e KS k £ C . (4.2) 

The equations for h p (x) are rather lengthy and can be found in [S] so we do not 
reproduce them here. Using the same strategy as in the previous section, we expand 
h p (x) and k in powers of k. The observed behavior is 



N p (x) = N p (x) , 
A p (x) = A p (x) , 



w p {x) 

!>l "' ; ~ ~~ w 



w p (x) 



v p (x) 



v p (x) 



upy ~' k 

K = R + 0{k 2 ) . 



(4.3) 
(4.4) 

(4.5) 

(4.6) 
(4.7) 



Critical Collapse of an Ultrarelativistic Fluid in the T — * 1 Limit 



8 



After substituting (|3.2U|) - (|3.23[1 and (|4.3|1 — (|4.T|I into the original equations and keeping 
only the leading terms we obtain a set of ODE's 
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The algebraic equation has the form 



(1 - n)A p - 2N e x uj(v v + v lo p ) 



2uouip = 



(4.8) 



(4.9) 



which can be used to verify the consistency of the results. 

The solution methods for the limiting case are the same as described in In 
particular to solve for the v(x) and u)(x), we must regularize the equations (|3.24|) - 
<|3.25fl at the sonic point (x = 0) and find the value of v(0) for which the solution has 
the correct behavior as x — > — oo. Similarly, the value of R is found by requiring that 
Vp{x) does not diverge as x — > — oo. To numerically integrate the set of ODE's we used 
the package LSODE \Y&\ I14j . For increased precision we used 128 bit representation 
of floating point numbers (quadruple precision). Figure [21 shows the limiting solutions 
for the relevant eigenmodes. 



5. Results of numerical calculations 



Numerical modeling of type II critical collapse of relativistic fluids even in spherical 
symmetry is a rather challenging task. The reason is that the dynamics takes place 
at ever decreasing spatial length scale and also the density and pressure of the fluid 
increase by many orders of magnitude. For example, in our calculations the density 
and pressure increase typically by a factor of million for the nearly critical solutions 
(in the subcritical regime) and in some of the supercritical calculations the central 
density reached a value greater than 10 54 . It is clear that some kind of adaptivity is 
required in order to perform these calculations. 

For r close to 1, i.e., in the regime we want to explore, another difficulty arises — 
tuning the critical parameter to about 15 digits (the limit for double precision floating 
point numbers) is not sufficient and calculations in quadruple precision are necessary. 
This prolongs the runtime of the code significantly because on "standard" workstations 
based on the x86 architecture quadruple precision is not natively supported by the 
processor. A more detailed description of some of the specifics of our numerical 
implementation can be found in section |S| 




Figure 2. Plot of the relevant eigenmodes of the limiting CSS solution N p (x), 
A p (x), ui p (x), v p (x). 



5.1. Comparison of numerical and CSS critical solutions 

If our numerical scheme works properly we should obtain the same critical solutions 
as we did by using the CSS ansatz (up to a truncation error). We must bear in 
mind, however, that the CSS critical solutions do not describe an asymptotically 
flat spacetime whereas the spacetime generated numerically is asymptotically 
Schwarzschild and therefore the solutions match only in a limited domain close to 
r = 0. In order to compare the solutions we have to translate the coordinate r used 
in numerical calculations into the self similar coordinate x in which the CSS solutions 
are cast. The relation is provided by the equation l|3.3fl . Since the relation of t„ to 
t is not known we use some distinct feature of the solution, e.g., local minimum or 
maximum to identify a particular r with a particular x. This allows us to calculate 

Figures show the comparison for k 2 = 10~ 6 . The numerical data were taken 
from the closest subcritical run just before the fluid dispersed. The numerical solutions 
agree very well with the CSS ones in a limited region as expected. The results for 
other values of k agree similarly well. 

Note that the lapse a used in the numerical calculations is not the same as the 
a from equation (|3.4(l . Therefore in order to compare the functions N(x) we must 
rescale one by a constant factor — we simply matched the leftmost data point from 
the numerical calculations with the corresponding one from the CSS solution. The 
dotted line shows the AMR hierarchy level. An increase in hierarchy level by one 
corresponds to a reduction of the cell size by factor of 2. 
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Figure 3. Comparison of CSS and numerical solution for A(x) for k 2 = 10 



5.2. Determining the scaling exponent 

For the numerical calculations we chose the initial density profile to be a Gaussian 

p = pexp[-(r-r c ) 2 /A 2 ] , (5.10) 

and the velocity v was set to zero (we used r c = and A = 0.2). The amplitude of 
the Gaussian p serves as the tunable parameter. Its critical value p* can be found by 
a bisection search. Once the critical parameter is obtained with sufficient precision it 
is possible to calculate the scaling exponent. 

We calculate the scaling exponent from subcritical runs. It has the obvious 
advantage that a black hole does not need to form in our spacetime (in our coordinate 
system it is not even possible). A further explanation will be given in the next section. 
The scaling relation for the trace of the stress-energy tensor has the form 

max(T%) = 3P- P Hp-p*r 2 ^ ■ (5-11) 
An accurate determination of the scaling exponent in practice is not as 
straightforward as it might seem. The basic idea is to perform a number of subcritical 
runs and then fit a straight line to the logarithm of the equation (|5.11|l . In our 
approach we do not use a predetermined value of p* but also optimize p* in order to 
get the best fit. In that sense the fit is non-linear. 

We would like to have our data points separated approximately evenly in the 
log \p — p* | coordinate but since we do not know the value of p* beforehand we perform 
a trial fit to obtain the value of p* and then use that value as a reference. 
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Figure 4. Comparison of CSS and numerical solution for log(7V(a;)) for k 2 = 
10" 6 . 

Another problem is to determine the interval of log \p — p*| over which we perform 
the fit. A very broad interval is likely to give inaccurate results since the scaling 
behavior is valid only in the vicinity of . On the other hand fitting over a very 
narrow interval extremely close to the might produce incorrect result as well. We 
use the following approach which we believe is sufficiently robust and provides a way 
to estimate the error as well. 

First we generate data points that span rather wide range of log \p— p* \ . Typically, 
we use around one hundred data points. Then we perform what we call a windowed fit, 
i.e., fit data points with the index i — k, . . . , k + N w — 1 where k = 1, . . . , N — N w + 1. 
N w is the width of the window and N is the total number of data points. If N w = N 
then the window spans all the data points. 

Figure shows data obtained from subcritical runs for k 2 = 10 -2 . The dashed 
line (not very visible) corresponds to the best fit and yields a value of Tat = 0.1148. 
This is the value we report. Figure |H1 shows 7 as a function of k for various window 
sizes (different markers correspond to different window sizes). Each data point is 
plotted with an error bar corresponding to the error of each 7 estimate that is equal 
to one half of the standard deviation corresponding to the slope The smaller 

the error bar the better the linear fit. The number of data points for each window 
size corresponds to the largest k of N — N w + 1. The main purpose of this plot is to 
provide some quantitative estimate of the error of the Tat . In addition, it may provide 
a hint whether there is some trend in 7 as we fit data closer to the critical value p*. 
Obviously, the largest variations are present for the smallest window sizes since we fit 
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Figure 5. Comparison of CSS and numerical solution for w{x) for k 2 = 10 
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Figure 6. Comparison of CSS and numerical solution for v(x) for k 2 = 10 6 . 
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9.463694624 
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Table 1. Scaling exponents calculated using both the CSS ansatz and a direct 
numerical solution. The values reported for k = are the limiting values. 



less data points. However, the averages 

^ N-N„ + l 

^= N-N W + 1 E 7k (5 ' 12) 

k— 1 

for different window sizes N w show very little variation. This is illustrated in figure^] 
Looking at figure |H1 (although rather messy) we can very conservatively declare 
that for k 2 = 10~ 2 

7 < 0.1153. (5.13) 

The plot also suggest a trend of decreasing 7 as we fit data closer to . Therefore it 
is not straightforward to estimate the lower bound for 7. 

We created the same type of plots for all the other values of k. Their structure is 
similar therefore we only show one additional set for k 2 = 10~ 5 ffigures HT!Hl2|) . For 
k 2 = 10 -5 we do not see any obvious decreasing trend for 7. Moreover, we observe a 
drop in the 7 value for the smallest window sizes plotted and large values of k. This 
just illustrates the fact that it could be dangerous to rely only on data from a narrow 
interval around the critical value . 

Table ^ summarizes the results obtained from both the perturbation theory (7 SS ) 
and the numerical calculations (7at)- The error reported is the percentage difference 
between 7 SS and "fat- Our results are in perfect agreement with [5] for the values of k 
reported therein. 

5.3. The supercritical regime 

In the previous section we showed how to calculate the critical exponents from 
subcritical runs using l|5.11|l . Originally, the order parameter was taken to be the 
mass of a black hole formed during the supercritical collapse. Why haven't we then 
used the supercritical runs and the black hole mass to determine the scaling exponent? 
Aside of the fact that in the coordinate system we are using we can not really form a 
black hole a more fundamental reason exists. There is no sign of black hole formation 
in the supercritical runs. Even after the central densities reach large values, there is no 
sign of event horizon formation (i.e. 2m /r approaches a constant value smaller than 
1). According to jHj this is expected and the universal attractor is not a spacetime with 
a black hole but a general relativistic Larson-Penston solution (GRLP). The GRLP 
solution exists only for T — 1 < 0.036 ± 0.002 and it contains a naked singularity for 

r - 1< 0.0105 [IT]. 
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Figure 7. Fitted data from subcritical solutions for k = 10 
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Figure 8. Windowed fits of data from subcritical runs for k 2 = 10 2 , AT, 
30,40,50,60,70,80,90. 
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Figure 9. Averaged 7 as a function of window size for k = 10" 
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Figure 10. Fitted data from subcritical solutions for k 2 = 10 5 . 
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Figure 11. Windowed fits of data from subcritical runs for k 2 = 10 
N w = 40, 50, 60, 70, 80, 90. 
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Figure 12. Averaged 7 as a function of window size for k 2 = 10 5 . 
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The GRLP is a solution of equations (|3.13[) - ([3.17[1 but unlike the critical solution 
(which is a general relativistic generalization of the Hunter's type (a) solution) it is a 
"pure collapse" solution, i.e., the velocity is always negative. 

To test this hypothesis, we performed a generic supercritical run for T — 1 = 0.01 
and r — 1 = 10~ 6 . For T — 1 = 10 -6 we stopped the calculations when the refinement 
level reached 100. At that point the central density reached value of about 10 54 and 
the spatial resolution at the center was about 10~ 32 . For T - 1 = 0.01 we stopped 
the calculations at refinement level 65. The central density reached value of about 
10 38 and the spatial resolution at the center was about 10~ 22 . We stress that the 
calculations were stopped artificially and they could be, if needed, pushed further. 
We also performed a control supercritical run for T — 1 = 0.02 and we observed the 
quantity 2m /r approaching 1 as expected. 

Figures ITTA and ITU show the comparison of the GRLP solution with the numerical 
data for r — 1 = 0.01. Figures IT51 and [TBI show the comparison of the GRLP solution 
with the numerical data for T — 1 = 10~ 6 . In the plots we show only a fraction of 
the data points outside of the very central region so the GRLP solution is visible on 
the plots. The agreement is very good near the origin and starts to deviate at larger 
distances. This is expected since the GRLP solution is not asymptotically flat whereas 
the numerical solution is. 

These results support the hypothesis that the GRLP solution is an attractor for 
the supercritical collapse and therefore naked singularities can be formed in a generic 
gravitational collapse. It is conceivable that black hole is also an attractor for the 
supercritical collapse. It would certainly be interesting to investigate the regime of 
the transition from GRLP to a black hole spacetime. 

6. Numerical scheme 

Our numerical scheme for the fluid evolution is based on high resolution shock 
capturing methods (HRSC) that have been used extensively in recent years. They are 
described in details in numerous papers and review articles [151 1161 1171 H%] . Therefore 
we are going to describe only those features of our code that are new or not completely 
obvious. All the calculations have been performed with quadruple precision arithmetic 
in order to tune the critical parameter p to as many as 30 significant digits. 

We used adaptive mesh refinement (AMR) in order to follow the evolution on an 
ever decreasing spatial and time scales. Another feature we implemented was a "true" 
vacuum, i.e., in the vacuum regions we do not keep a small residual fluid atmosphere 
(a so called floor) but the fluid variables are set to zero values. 

The regularity conditions for fluid variables are achieved by using reflective 
boundary conditions at the origin r = 0. We use outflow boundary conditions at 
the outer boundary (although there should be no fluid crossing the outer boundary). 

To solve for the geometry we first solve for a using the equation l[2.27|) . We 
demand that a{r = 0) = 1 and integrate the equation outward, a is obtained from 
equation l[2.26[) by integrating it inward. It is initialized by setting a(r = r max ) = 
-l/a(r = r max ). 
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Figure 13. Comparison of the supercritical numerical solution and the GRLP 
solution for T - 1 = 0.01. 

6.1. Implementation of AMR 

Our implementation of AMR can be viewed as an evolution on a dynamically changing 
non- uniform grid with hierarchical structure!. We use a single global time step that 
is determined by the size of the smallest cell. 

The grid structure changes if certain criteria are not met. The tests are performed 
at every fixed number of steps (at a checkpoint). The criteria we first implemented 
were based on the truncation error estimate. We compare a one step update on a 
2:1 coarsened grid with two step update on the original grid. The problem with this 
approach was that occasionally we experienced problems with the update during the 
coarse grid step. This certainly could be fixed but we opted for a simpler criteria. In 
essence we want to guarantee that all the features of the solutions are well resolved. 
We demand higher resolution in regions with larger differences in gradients of the 
variables. The cells that do not meet our criteria are flagged (we can use more than 
one criterion). 

After flagging we apply buffering, i.e., we flag nearby regions around the originally 
flagged cells. The size of the buffer regions is chosen so that each of the flagged cells 
could affect at most the cells within the buffer region during the evolution until the next 
checkpoint. Or, equivalently, we can say that only information from the buffer region 
could have reached the flagged cell during evolution from the previous checkpoint. 

After the buffering step the grid is adjusted so that all the flagged cells are refined 
(we use 2:1 refinement ratio). We demand that the refinement algorithm preserve the 

\ In this context a hierarchical grid is a grid in which the size of two neighboring cells may differ at 
most by the refinement ratio (in our case 2). 
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hierarchical structure. 

At this point we restart the evolution from the last successful checkpoint but with 
the newly created grid. All the variables are interpolated onto the new grid. If no 
cells were flagged (the grid is not changed) then the current state is stored (the grid 
structure and all the variables) and the evolution continues. 

In systems where discontinuities may be present we must be a little bit more 
careful with the refinement scheme. Typically, our criteria would always fail at 
discontinuities and therefore a naive application of the aforementioned rules would 
lead to an uncontrolled refinement (in practice, we always set a maximum refinement 
level). This is not desirable since this would slow down the evolution tremendously. 
Moreover, our scheme is "shock capturing" therefore it should treat discontinuous 
solutions properly. In general this issue is not trivial an some clever scheme must be 
applied to deal with this. Our situation is slightly easier since we qualitatively know 
the dynamics. Although we do not expect shocks in the critical solutions we do have 
shocks in the computational domain. These are present in the "outer domain", i.e., 
beyond the central self similar region. We therefore employ a simple approach - 
we introduce different refinement maxima for the inner and outer regions (the outer 
region maximum is much smaller than the inner one). 

6. 2. Treatment of the vacuum 

Vacuum regions always pose a problem in numerical hydrodynamics. Technically, the 
problem is that in very rarefied regions one part of the numerical scheme fails — in 
particular the conversion from conservative to primitive variables. It simply happens 
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Figure 15. Comparison of the supercritical numerical solution and the GRLP 
solution for T - 1 = 10~ 6 . 

that there are no physical values of the primitive variables that correspond to the 
updated conservative variables. These regions are, however, dynamically completely 
unimportant since the energy and momentum densities there are negligible. 

A standard approach is to maintain a minimal "atmosphere" of fluid everywhere 
with the density several orders of magnitudes smaller than the typical densities in the 
system. Although we believe that there is nothing wrong with this approach and it 
would most certainly work in our situation we tried to develop a different approach. 

One nice feature of a "true" vacuum is that there are no fluxes through the outer 
boundary and therefore no artificial reflections. 

Of course we can not really solve the problem because of the non-existence of the 
exact Ricmann solution at the vacuum boundary region. What we propose is a set of 
ideas which typically work in practice. They are not universal and one might need to 
adjust them slightly for different systems. Sometimes not all of them need to be used. 
In general, the situation is more complicated with stiffer fluids (T — > 2) but we were 
able to use the scheme even for critical collapse of an ideal fluid with T = 2. We also 
used this approach in 2D modeling of axisymmetric fluid accretion onto black holes 
and some preliminary 3D runs of neutron stars. For the type of systems studied in 
this paper (r — > 1) the scheme works in its simplest form. 

The basic ideas are as follows. 

(i) set the flux between two vacuum cells to zero 

(ii) the flux between vacuum and non-vacuum cell is calculated using the Tadmor's 
scheme (as described in 




Figure 16. Comparison of the supercritical numerical solution and the GRLP 
solution for T — 1 = 10~~ 6 — detail of the central region. 



(iii) after the update of the conservative variables set to vacuum all cells with fluid 
levels below certain threshold 

The purpose of item (|IjJ is to use a simple and robust update scheme which relies 
only on the characteristic speeds and not on the full spectral decomposition typically 
needed in Roe type of solvers. 

Often, the approach described above is sufficient. This was the case for all the 
calculations we performed for purpose of this paper. Since we used quadruple precision 
and our fluid is extremely soft we did not even need to apply the rule number 2. 

Sometimes we run into problems even after applying the above ideas. In this case 
we should always try to adjust the vacuum threshold up to the maximum level we can 
afford without loosing a significant amount of fluid in the system. 

If we use AMR and we run into trouble at some cell we can flag the cell and 
invoke an "emergency" refinement and restart, especially if the cell is not at the 
vacuum boundary. We used this approach when we tested the numerical scheme on 
r = 2 ideal fluid. If the difficulty arises at the vacuum boundary (which happens in 
the majority of cases) we can simply set the cell to the vacuum values. 

In general, it is likely that the velocities near the vacuum boundary will not be 
smooth and will be be highly relativistic but this should not have any adverse effect 
on the dynamics (as long as the scheme works). 
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7. Conclusions 



In this paper we calculated the critical solutions and scaling exponents for values of 
T close to 1 using both the CSS ansatz and by numerically solving the full set of 
equations for spherically symmetric ultrarelativistic fluid collapse. Moreover, by using 
the CSS ansatz we obtained and solved the equations for the limiting case T = 1. In 
the limit of T — ► 1 the general relativistic critical solution converges to the Newtonian 
Hunters type (a) solution. The limiting value of the scaling exponent 

lim 7(» = 0.1056669768 (7.14) 

can be therefore interpreted as the scaling exponent for critical collapse in Newtonian 
gravity. This is consistent with conclusions drawn from previous work |1 II [U] . 

We also addressed the possibility of generic naked singularity formation in the 
supercritical regime. Our calculations support the fact that the supercritical collapsing 
solutions for T—l < 0.01 converge to the GRLP solution (an in the limit of T — > 1 to the 
Newtonian Larson-Penston solution) and from then follows that these solutions 
contain a naked singularity. In other words, the GRLP solution is an endstate for 
a supercritical collapse for a slightly supercritical dataset for our family of initial 
data§. It would be interesting to further explore the structure of the solution space, 
in particular, to search for a supercritical collapse with black hole as an endstate. 

The numerically found solutions and scaling exponents agree very well (to a 
fraction of a percent) with those obtained by using the CSS ansatz. To numerically 
solve the equations we used advanced numerical methods including AMR, special 
vacuum treatment and quadruple precision floating point arithmetic. 
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